---
jupytext:
  formats: md:myst,ipynb
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.4
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

```{code-cell}
:init_cell: true
:tags: [hide-cell]

import mmf_setup; mmf_setup.nbinit()
```

# Exact Solutions

Here we present the various exact solutions described in
[`gpe/exact_solutions.py`](../gpe/exact_solutions.py)

## Harmonic Oscillator

+++

Here we construct an exact solution for a trapped BEC in a harmonic oscillator:

\begin{gather*}
  \psi_0(x) = \sqrt{n_0}e^{-x^2/2\sigma^2}, \qquad
  V(x) = \frac{m\omega^2x^2}{2} - \frac{\hbar^2}{2m\sigma^2} - gn_0e^{-x^2/\sigma^2}, \qquad
  \sigma = \sqrt{\frac{\hbar}{m\omega}}.
\end{gather*}

This potential has the non-linear piece subtracted so that the constructed solution is an exact solution with zero chemical potential (i.e. a stationary state) and energy:

\begin{gather*}
  E = \int_{-\infty}^{\infty}\d{x}\;\frac{-gn_0^2}{2}e^{-2x^2/\sigma^2} = -\sqrt{\frac{\pi}{2}}\frac{gn_0^2\sigma}{2}
\end{gather*}

```{code-cell}
%pylab inline --no-import-all
from gpe.imports import *
from gpe.minimize import MinimizeState
from gpe.exact_solutions import HarmonicOscillator
```

```{code-cell}
s = HarmonicOscillator(n0=1.0)
s.plot()
assert abs(s.compute_dy_dt(s.copy())[...]).max() < 1e-14
assert np.allclose(s.get_energy(), -s.g * s.n0**2 * s.sigma * np.sqrt(np.pi / 2) / 2)
```

We can also do this in higher dimensions:

\begin{gather*}
  \psi_0(\vect{x}) = \sqrt{n_0}\exp\left(-\sum_{i}\frac{x_i^2}{2\sigma_i^2}\right), \qquad
  \frac{-\hbar^2}{2m}\nabla^2\psi_0 
  = \frac{\hbar^2}{2m}\sum_{i} \left(\frac{1}{\sigma_i^2} - \frac{x_i^2}{\sigma_i^4}\right)
  = \sum_{i} \left(\frac{\hbar^2}{2m\sigma_i^2} - \frac{m\omega_i^2x_i^2}{2}\right)
  \\
  \qquad
  V(x) = \sum_i\left(\frac{m\omega_i^2x_i^2}{2} - \frac{\hbar^2}{2m\sigma_i^2}\right) - g\abs{\psi_0}^2, \qquad
  \sigma_{i}^2 = \frac{\hbar}{m\omega_i}.
\end{gather*}

This potential has the non-linear piece subtracted so that the constructed solution is an exact solution with zero chemical potential (i.e. a stationary state) and energy:

\begin{gather*}
  E = \int_{-\infty}^{\infty}\d{x}\;\frac{-gn_0^2}{2}e^{-2x^2/\sigma^2} = -\sqrt{\frac{\pi}{2}}\frac{gn_0^2\sigma}{2}
\end{gather*}

+++

# Travelling Waves

+++

The standard GPE admits a family of periodic traveling waves with an analytic form:

\begin{gather*}
  \psi(x, t) = e^{-\I \mu t /\hbar}\psi_v(x_v), \qquad
  x_v = x - vt
\end{gather*}

where the wave is moving with speed $v$.  We base our solution on the presentation of Hoefer and El [El:2016].

The GPE can be expressed in terms of $\psi(x, t)$ or in terms of $\psi_v(x_v) = e^{-\I\mu t/\hbar} \psi_v(x- vt)$:

\begin{gather*}
  \I\hbar\dot{\psi} = -\frac{\hbar^2}{2m}\pdiff[2]{\psi}{x} + g\abs{\psi}^2\psi,\qquad
  0 = -\frac{\hbar^2\psi_v''}{2m} + v \I \hbar \psi_v' + (g\abs{\psi_v}^2 - \mu)\psi_v
    = \left(\frac{\op{p}^2}{2m} - v\op{p} + g\abs{\psi_v}^2 - \mu\right)\psi_v.
\end{gather*}

*Note: this is not the typical Galilean transformation for quantum mechanical systems which includes an extra factor of $e^{\I (mvx +mv^2t/2)\hbar}$.  This additional phase allows one to shift the momentum, completing the square of the kinetic energy, and removing the $v\op{p}$ linear term at the expense of shifting the chemical potential.  The present form, however, is more general, and works with arbitrary dispersion, so we maintain it.*

Hoefer and El [El:2016] set $\hbar=m=g=1$ which defines the following dimensions:

\begin{gather*}
  m_\mathrm{unit} = m, \qquad
  l_\mathrm{unit} = \frac{mg}{\hbar^2}, \qquad
  t_\mathrm{unit} = \frac{m^3g^2}{\hbar^5}.
\end{gather*}

Using the Madelung transform:

\begin{gather}
  \psi(x, t) = \sqrt{\rho(x, t)}e^{\I\phi(x, t)}, \qquad u(x, t) = \phi'(x, t),\\
  \dot{\rho} + (\rho u)' = 0, \qquad
  (\rho u)_{,t} + \left(\rho u^2 + \frac{\rho^2}{2}\right)'
  = \frac{1}{4}\bigl(\rho\, (\log\rho)''\bigr)'. \tag{2.111}
\end{gather}

The solution can be expressed as:

\begin{gather*}
  \DeclareMathOperator{\sn}{sn}
  \psi(x, t) = e^{-\I\mu t/\hbar}\psi_v(x_v), \qquad
  \psi_v(x_v) = \sqrt{\rho_v(x_v)}e^{\I\phi_v(x_v)}, \qquad
  u(x, t) = \phi_v'(x, t),\\
  m = al^2, \qquad
  \rho_v(x_v) = \rho_{\min} + a \sn^2\bigl(\frac{x_v}{l};m\bigr), \qquad
  x_v = x - vt, \qquad
  u(x, t) = v - \frac{C}{\rho(x, t)}, \\
  a = \rho_{\max} - \rho_{\min}, \qquad
  \mu = \rho_{\min} - \frac{v^2}{2} +  \frac{\rho_\max + a/m}{2}\qquad
  C = \frac{\sqrt{\rho_{\min}\rho_{\max}(1+\rho_{\min}l^2)}}{l},\qquad
  L = 2lK(m), \qquad
  Q = \frac{2\pi}{L}.
\end{gather*}

*(Note: $k=\sqrt{m}$ for $\sn(z,k)$ in some CASs.  We use $\sn(z;m)$ and $K(m)$ here)*

[El:2016]: https://doi.org/10.1016/j.physd.2016.04.006 'G. A. El and M. A. Hoefer, "Dispersive shock waves and modulation theory", Physica D: Nonlinear Phenomena 333, 11--65 (2016) [arXiv:1602.06163](http://arXiv.org/abs/arXiv:1602.06163)'

The full solution (with proper coefficients) is thus:


\begin{gather*}
  \psi(x, t) = e^{-\I\mu t/\hbar}\psi_v(x_v), \qquad
  \psi_v(x_v) = \sqrt{\rho_v(x_v)}e^{\I\phi_v(x_v)}, \qquad
  u(x_v) = \frac{\hbar}{m}\phi_v'(x_v) = v - \frac{C}{\rho_v(x_v)},\\
  \rho_v(x_v) = \rho_{\min} + a\sn^2\bigl(\frac{x_v}{l};\tfrac{mg}{\hbar^2}al^2\bigr), \qquad
  x_v = x - vt, \qquad
  a = \rho_{\max} - \rho_{\min}, \\
  \mu = g\rho_\min - \frac{mv^2}{2} + g\frac{\rho_\max +a/m}{2}\qquad
  C = \frac{\sqrt{\rho_{\min}\rho_{\max}(\hbar^2+mgl^2\rho_{\min})}}{ml}.
\end{gather*}

+++

\begin{gather*}
  m = al^2, \qquad
  l^{-1} = \frac{2K(m)}{L}
\end{gather*}

```{code-cell}
%pylab inline --no-import-all
from gpe.imports import *
import gpe.exact_solutions

reload(gpe.exact_solutions)
from gpe.exact_solutions import TravellingWaves

# Stationary wave
args = dict(n0=0.9, n1=1.0, Lx=10.0, Nx=16)
s0 = TravellingWaves(v_p=0, **args)

# Travelling wave in the lab frame
s1 = TravellingWaves(v_p=0.6047197603, v_x=0, **args)

# Travelling wave in a co-moving frame.
s2 = TravellingWaves(v_p=0.6047197603, **args)

for ss in evolves([s0, s1, s2], t_max=50.0):
    plt.clf()
    for s in ss:
        s.plot()

# s.twist_x
# s.mu, s.get_mu().real, s._twist
```

As a simple test, we consider the phonon limit $a = \epsilon \ll 1$.  In this limit, we have

\begin{gather*}
  K(m) \approx \frac{\pi}{2}\left(1 + \frac{m}{4}\right) + \order(m^2), \qquad
  m = \frac{L^2}{\pi^2}\epsilon, \qquad
  \sn(\theta, m) = \sin(\theta) + \order(m), \qquad
  l = \frac{L}{\pi}\left(1 - \frac{L^2}{4\pi^2}\epsilon\right) + \order(\epsilon^2)\\
\end{gather*}

\begin{gather*}
  \frac{L = \approx l\pi\left(1 + \frac{m}{4}\right)
\end{gather*}

consider a $v=0$ solution with $\rho_\min = \rho_\max = 1$:

\begin{gather*}
  \rho = 1, \qquad a = m = 0, \qquad
  \frac{L}{l} = \pi, \qquad
  u(x_v) = \frac{C}{\rho}, \qquad
  \theta(x_v) = \frac{C}{\rho}x_v
\end{gather*}

```{code-cell}
from gpe.exact_solutions import K
```

```{code-cell}
s._C
```

```{code-cell}
s = TravellingWaves(n0=1.0, n1=1.0, Lx=np.pi, Nx=64, v_p=0.0)
s.plot()
s.twist_x
s._mu, s.get_mu().real, s._twist
```

```{code-cell}
e = EvolverABM(s, dt=0.2 * s.t_scale)
with NoInterrupt() as interrupted:
    while not interrupted:
        e.evolve(1000)
        plt.clf()
        e.y.plot()
        display(plt.gcf())
        plt.close("all")
        clear_output(wait=True)
```

```{code-cell}
x = s.xyz[0]
c = np.sqrt(s.g * s.n1 / s.m)
v = 0.5 * c
u = np.sqrt(c**2 - v**2)
n0 = s.n1 * (v**2 / c**2)
l = s.hbar / s.m / u
psi_soliton = v / c * 1j + u / c * np.tanh(x / l)
twist = np.angle(psi_soliton)[-1] - np.angle(psi_soliton)[0]
s1 = TravellingWaves(n0=n0, n1=s.n1, Lx=s.basis.Lx, Nx=s.basis.Nx, v_p=0, twist=twist)
s1[...] = psi_soliton
s1.plot()
```

```{code-cell}
e = EvolverABM(s, dt=0.2 * s.t_scale)
e1 = EvolverABM(s1, dt=0.2 * s1.t_scale)
pe = None
with NoInterrupt() as interrupted:
    while not interrupted:
        e.evolve(1000)
        e1.evolve(1000)
        plt.clf()
        e.y.plot()
        e1.y.plot()
        display(plt.gcf())
        plt.close("all")
        clear_output(wait=True)
```

(sec:BrightSoliton)=
# Bright Soliton

+++

Here we demonstrate the analytic bright-soliton for a GPE with attractive interactions.  This object moves at a specified speed $v$, and we consider the solution in three frames: comoving (the soliton is stationary), lab frame (the soliton makes a single oscillation in time $T=L/v$) and a frame moving backwards at the same speed (the soliton crosses the box twice in time $T$).

*This forms the basis of [`gpe/tests/test_bec.py::test_StateTwist_x_v_x`](../gpe/tests/test_bec.py) and verifies that the frame velocity argument of `StateTwist_x` works.  Because the density vanishes at the boundaries, this does not test the twist.*

+++

\begin{gather*}
  \I\hbar\dot{\psi}(x, t) = \frac{-\hbar^2\psi''(x, t)}{2m} + g\abs{\psi}^2\psi,\\
  \psi(x, t) = \frac{\sqrt{n_0}}{\cosh\bigl(\eta (x-vt)\bigr)}
  \exp\left\{\frac{\left(\mu + \tfrac{mv^2}{2}\right)t - vx}{\I\hbar}\right\}, \qquad
  \mu = \frac{-\hbar^2\eta^2}{2m}, \\
  gn_0 = 2\mu = \frac{-\hbar^2\eta^2}{m}, \qquad
  \sigma = \frac{1}{\eta}.
\end{gather*}

```{code-cell}
import sympy
from sympy import Eq, var, I, exp, cosh, sqrt, Abs
```

```{code-cell}
x, t, g = var(["x", "t", "g"], real=True)
eta, n0, mu, m, hbar = var(["eta", "n_0", "mu", "m", "hbar"], positive=True)
mu = -(hbar**2) * eta**2 / 2 / m
g = 2 * mu / n0
psi = exp(mu * t / I / hbar) * sqrt(n_0) / cosh(eta * x)
n = Abs(psi) ** 2
display(Eq(sympy.S("psi"), psi))
Hpsi = -(hbar**2) * psi.diff(x, x) / 2 / m + g * n * psi
display(Eq(sympy.S("H*psi"), Hpsi))
(Hpsi - I * hbar * psi.diff(t)).simplify()
```

```{code-cell}
%pylab inline --no-import-all
from gpe.imports import *
import gpe.exact_solutions

reload(gpe.exact_solutions)
from gpe.exact_solutions import BrightSoliton

Nx = 32
Lx = 10.0
v = 2.0

args = dict(Nx=Nx, Lx=Lx, v=v, sigma=1)
s0 = BrightSoliton(v_x=0, **args)
s1 = BrightSoliton(v_x=-v, **args)
s2 = BrightSoliton(v_x=v, **args)

for ss in evolves([s0, s1, s2], t_max=Lx / v):
    plt.clf()
    for s in ss:
        s.plot()
```

## 1D GPE (NLSEQ)

:::{margin}
I do not understand this yet.  It is also called a Zakharov-Shabat system.
:::
Here we consider the exact solutions to the 1D GPE via the inverse scattering method
(ISM). These follow from a so-called Lax representation, which can be expressed as
follows.  Consider two linear operators $\op{L}(\lambda)$ and
\begin{gather*}
  \op{M}(\lambda) = \I\diff{}{t} + V(x, t, \lambda)
\end{gather*}
that commute:
\begin{gather*}
  [\op{L}(\lambda), \op{M}(\lambda)] = 0.
\end{gather*}

\begin{gather*}
  \mat{L} = \I\mat{1}\pdiff{}{x} + \begin{pmatrix}
    0 & u^*\\
    u & 0
  \end{pmatrix},\\
  \mat{M} = -\mat{1}\pdiff[2]{}{x} + \begin{pmatrix}
    \abs{u}^2 & \I u_{,x}^*\\
    -\I u_{,x} & -\abs{u}^2
  \end{pmatrix}.
\end{gather*}









differential equation
\begin{gather*}
  \I \dot{\mat{L}} = [\mat{L}, \mat{A}]
\end{gather*}



:::{margin}
We can obtain this form by choosing units so that
\begin{gather*}
  2m = \hbar = 1,
\end{gather*}
and scaling
\begin{gather*}
  u^2 = \frac{2mg}{\hbar^2}\psi^2.
\end{gather*}
:::
We start by expressing the problem as
\begin{gather*}
  \I \dot{u} + u'' + 2 \abs{u}^2 u = 0, \qquad
  u(x, t=0) = u_0(x).
\end{gather*}
Define
\begin{gather*}
  \mat{q}_0(x) = \begin{pmatrix}
    0 & u_0(x)\\
    u_0^*(x) & 0
  \end{pmatrix}.
\end{gather*}